Flowing sand - a physical realization of Directed Percolation 
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We introduce and investigate a simple model to describe recent experiments by Douady and Daerr 
on flowing sand. The model reproduces experimentally observed compact avalanches, whose open- 
ing angle decreases linearly as a threshold is approached. On large scales the model exhibits a 
crossover from compact directed percolation to directed percolation; we predict similar behavior for 
the experimental system. We estimate the regime where "true" directed percolation morphology 
and exponents will be observed, providing the first experimental realization for this class of models. 



PACS numbers; 45.70.Ht, 64.60. Ht, 64.60.Ak 

Directed Percolation (DP) is perhaps the simplest 
model that exhibits a non-equilibrium phase transition 
between an "active" or "wet" phase and an inactive "dry" 
one . In the latter phase the system is in a single "ab- 
sorbing" state; once it reaches the completely dry state, 
it will always stay there. 

Interest in DP stems mainly from universality of the as- 
sociated critical behavior. It is believed that transitions 
in all models with an absorbing state are in the DP uni- 
versality class (unless there are some special underlying 
symmetries). Even though DP exponents have not yet 
been calculated analytically, their values were measured 
for a wide variety of models and are known (especially in 
1-1-1 dimensions) with very high precision B. 

Models in the DP universality class are supposed to de- 
scribe a wide variety of phenomena, ranging from catal- 
ysis to turbulence |^; nevertheless, so far no physical 
system has been found that exhibits DP behavior and 
exponents We show here that a simple system of 
sand flow on an inclined plane, recently studied by Daerr 
and Douady (DD) |], may well be the first physical re- 
alization of a transition in the DP universality class. 

The results reported by DD present a puzzle, namely 
the threshold phenomenon they discovered exhibits 
"wet" clusters whose shapes differ from those seen in 
standard DP simulations; they are much more compact. 
The corresponding model, called Compact Directed Per- 
colation (CDP) is unstable against perturbations towards 
the standard DP behavior Q; the latter is the generic 
case expected to occur. Since DD did no fine-tuning to 
place their system in the CDP class, their observations 
are surprising. 

This motivated us to look for a simple model, which 
is defined in terms of dynamic rules that can be plausi- 
bly related to the experiments and exhibits features that 
reproduce the experimentally observed ones. We then 
investigated whether the transition exhibited by such a 
model belongs to the DP universality class. 

We propose a directed sandpile model, simpler than 
that of Tadic and Dhar Q ; here after each avalanche the 
system is reset to a uniform initial state. Our model has 



a transition from an inactive to an active phase, in which 
we see avalanches whose compact shapes reproduce the 
experimental observations and indeed, the resulting crit- 
ical behavior is close to CDP, rather than to DP. We re- 
solve this by showing that the CDP type critical behavior 
is a transient: the true critical behavior is of the DP type, 
but it can be seen only after a very long crossover regime. 
Our conclusion is that the DD experiment does serve as a 
possible realization of a DP-type transition. We propose 
here ways to shorten the crossover regime and to extend 
the scale on which the experiments are performed. 

1. The Douady-Daerr experiment. Glass beads 
("sand") of diameter 250-425 /im [|| are poured uni- 
formly at the top of an inclined plane (size Iti), cov- 
ered by a rough velvet cloth; the angle of inclination ipQ 
can be varied. As the beads flow down, a thin layer of 
thickness h = hdiifo), consisting of several monolayers, 
settles and remains immobile. At this thickness the sand 
is dynamically stable; the value of decreases with in- 
creasing angle of inclination. 

For each ipo there exists another thickness with 
hs{ipo) > hd{ipo), beyond which a static layer becomes 
unstable. Hence there exists a region in the (tp, h) plane, 
in which a static layer is stable but a flowing one is unsta- 
ble. We can now take the system, that settled at hd(}po), 
and increase its angle of inclination to ip = (po + Aip, stay- 
ing within this region of mixed stability. The layer will 
not flow spontaneously, but if we disturb it, generating a 
flow at the top, an avalanche will propagate, leaving be- 
hind a layer of thickness hd{(p). These avalanches had the 
shape of a fairly regular triangle with opening angle 6. 
As the increment Acp decreases, the value of 9 decreases 
as well, vanishing as Aip 0. This calls for testing a 
power law behavior of the form 



(1) 



If instead of increasing tp we lower the plane, i.e., go to 
Ap < 0, the thickness of our system, hdi^po), is less than 
the present thickness of dynamic stability, hd{^p). In this 
case an initial perturbation should not propagate, it will 
rather die out after a certain time (or beyond a certain 
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FIG. 1. The top layer of sand in a washboard potentiaL 

size of the transient avalanche). As \Aip\ 0, we 
expect this decay length to grow with a power law: 

^ i^A^r-^ii . (2) 

Hence by pouring sand at inclination ipQ, DD produced a 
critical system, precisely at the borderline (with respect 
to changing the angle) between a stable regime ip < (pQ, in 
which perturbations die out, and an unstable one (p > tpo, 
where perturbations persist and spread. The preparation 
procedure can be considered as a special kind of self- 
organized criticality (SOC) which differs from standard 
SOC models js) in which a slow driving force (acting 
on a time scale much smaller than that of the system's 
dynamic response) causes evolution to a critical state. 
Here avalanches are started by hand one by one. 

To associate this threshold phenomenon with DP, de- 
note by p the percolation probability and by pc its critical 
value. We associate the change in tilt with p — pc, i.e., 
assume that near the angle of preparation lpq the behav- 
ior of the sand system is related to a DP problem with 
A(p cx p — Pc- The exponent should be compared with 
the known values for DP and CDP. The exponent x in 
Eq. (|l|) can also be measured and compared with 

tanO^U/^W-iA^r^-'^ . (3) 

2. Definition of the model. To write down a simple 
model based on the physics of the flowing sand, we adopt 
an observation made by DD, that in the regime of inter- 
est ((/? « ipo) grains of the top layer rest on grains of 
the layers below (rather than on other grains of the top 
layer). Hence the lower layers provide for the top one a 
washboard potential, as shown in Fig. 0. 

We place the grains of the top layer on the sites of a 
regular square lattice with row index t and columns i. At 
any given time a grain G may become active if at least 
one of its neighbors from the row above has been active 
at the previous time step. If AE{G), the total energy 
transferred from these neighbors, exceeds the barrier Ei, 
of the washboard, G becomes active, "rolls down" and 
collides with the grains of the next row. The energy it 
brings to these collisions is 1 -t- AE{G), where 1 is our 
unit of energy, representing the potential energy due to 
the height difference between two consecutive rows (see 
Fig. |l|). A fraction / of its total energy is dissipated; the 
rest is divided stochastically among its three neighbors 
from the lower row. 



The model is defined in terms of two variables; an ac- 
tivation variable Sj = 0, 1 and an energy Ej. The index t 
denotes rows of our square lattice and time; at time t we 
update the states of the grains that belong to row t. The 
model is controlled by two parameters: Eb, the barrier 
height, and /, the fraction of dissipated energy. 

The dynamic rules of our model are as follows. For 
given activities 5** and energies Ej we first calculate the 
energy transferred to the grains of the next row t + 1. To 
this end we generate for each active site three random 
numbers, zl{6) (with 6 = ±1,0) that add up to 1. The 
energy transferred to grain (t + 1, i), given by 

^El+' = (1-/) E SlsElszlsiS), (4) 

5=±1,0 

determines its activation: 

f+i ^ r 1 active if AEl+^ > Et , , , 

\ inactive if AEI+^ < Et . ^ ' 

Then the energies of the next row of grains are set: 

El+' = Sl+' il + AEl+'). (6) 

The three random numbers zl{S) represent the fraction 
of energy transferred from the grain at site (t, i) to the 
one at {t + l,i + 5). We add up the energy contributions 
from these active sites; the fraction 1 — / is noi dissipated; 
if the acquired energy Ai?*"'"^ exceeds E^, site {t + l,i) 
becomes active, rolls over the barrier and brings to the 
collisions (at time t -\-2) the acquired energy calculated 
above and its excess potential energy (of value 1). 

3. Qualitative discussion of the transition and con- 
nection to the experiment. Let us vary E^ at a fixed 
value of the dissipation. For small values of Ei, an ac- 
tive grain will activate the grains below it with high 
probability; avalanches will propagate downhill and also 
spread sideways. For a strongly localized initial activa- 
tion we should, therefore, see triangular shaped activated 
regions. As Ejj increases, the rate of activation decreases 
and the opening angle of these triangles should de- 
crease, until Eb reaches a critical value E^, beyond which 
initial activations die out in a finite number of time steps 
(or rows). These expectations are indeed borne out by 
simulations of the model: the dependence of E^ on the 
dissipation / is shown in Fig. |^. 

The physics of the process is captured by a simple 
mean-field type approximation, in which all stochastic 
variables are replaced by their average values. Consider 
an edge separating an active region from an inactive one. 
At time t sites to the left of i and i itself are wet, whereas 
i -\- l^i -\- 2, ... are dry. Will site i be wet or dry at the 
next time step? In our mean-field estimate of the an- 
swer, assuming that all wet sites at time t have the same 
energy i?*, the energy delivered to site i at time t -\- 1 
is A£:*+i = |(1 - /)(! + AS*), where we set in Eq. (|) 
all z{5) = 1/3. At the critical point we expect all en- 
ergies to just suffice to go over the barrier; hence set 
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FIG. 2. Phase diagram of the model for flowing sand. The 
full line represents the phase transition line. The mean field 
approximation of Eq. (j^ is shown as a dotted line. The inset 
shows the opening angle tan 9 as a function of — Efj. 

AE*^'^ = AE* — E^. Solving the resulting equation 
yields a rough estimate of the transition line, 

E^, = 2(1 - /)/(! + 2/) , (7) 

as shown in Fig. It is easy to produce better mean- 
field type estimates of the transition and to compute the 
corresponding energy profile of the wet region |^ . 

To connect our model to the DD experiments note 
that the tilt angle (p tunes the ratio between the barrier 
height and the difference of potential energies between 
two rows. When the system is prepared at ipo, this ra- 
tio is precisely E^^. When one increases the tilt angle 
to if > ipQ, El, (measured in units of the potential differ- 
ence) decreases and we have Ef, < E^. As the tilt angle is 
now reduced, the size of Ef, increases, until it reaches its 
critical value precisely at (fQ. Thus increasing E'f, in the 
model corresponds to lowering the tilt angle towards (/Sg 
where the system is precisely at its boundary of dynamic 
stability. 

Hence to reproduce the experiment we were looking 
for (a) fairly compact triangular regions of activation for 
El, < E^, and (b) an opening angle of these triangles 
which should go to zero as Ef, approaches E^ from below. 

We simulated the model defined in Eqs. (^-(^ and 
found that it indeed reproduces these qualitative features 
of the experiment (see Fig. ||). The avalanches shown 
were produced for dissipation / — 0.5, activating a single 
site at t = 0, to which an initial energy of Eq = 100 
was assigned. As long as Et, was not too close to E^ the 
observed avalanches were compact, triangular and with 
fairly straight edges. The edges became rough only very 
close to E^^, such as the one shown on the right hand 
side of Fig. ^. The opening angle of the active regions 6 
decreased as E^ increased towards E^, as indicated in in 
the inset of Fig. 0. From these simulations we estimated 




FIG. 3. Typical avalanches starting from a single seed 
with dissipation / — 0.5 far away and close to criticality. 



E^ and the exponent (see Eq. (||)) 

X = - i/_L = 0.98(5) ~ 1 . (8) 

The linear variation of tan(0) with Aip is in agreement 
with experimental measurements Our findings have 
to be compared with the mean-field theory suggested in 
Ref. [0 which predicts a square root behavior. 

4- Crossover to directed percolation. The linear law 
(^ is consistent with the critical exponents of CDF |l^ 

= 2 , i^±^l, P^O. (9) 

These observations pose, however, a puzzle: since one be- 
lieves that DP is the generic situation, one would expect 
to find non-compact active regions and DP exponents. 
In fact, according to the DP conjecture jl^ any continu- 
ous spreading transition from a fluctuating active phase 
into a single frozen state should belong to the univer- 
sality class of directed percolation (DP), provided that 
the model is defined by short range interactions with- 
out exceptional properties such as higher symmetries or 
quenched randomness. The present model has neither 
special symmetries nor randomness JTsI ; it has a fluctu- 
ating active phase and exhibits a transition, character- 
ized by a positive one-component order parameter, into 
a single absorbing state. Hence the phase transition of 
our model should belong to the DP universality class. 

In order to understand this apparent paradox we per- 
formed high-precision Monte-Carlo simulations for dissi- 
pation / = 0.5 (see 1^ for further details). We performed 
time-dependent simulations p^ , i.e., we toppled a sin- 
gle grain in the center of the top row and measured the 
survival probability P{t) and the number of active sites 
N{t). At criticality, these quantities exhibit an asymp- 
totic power law behavior 

P{t) - t~^ , N{t) - t" . (10) 

In the case of CDP these exponents are given by [|6|jll|] 
6 = 1/2 and = 0' whereas DP is characterized by 
the exponents |^ J = 0.1595 and ry = 0.3137. Detecting 
deviations from power-law behavior in the long-time limit 
we estimated E^ = 0.385997(5). 
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FIG. 4. Average number of active sites N(t) ^^id mean 
survival probability P(t) of the toppling process at criticality 
averaged over 50 000 runs. The predicted slopes for CDP and 
DP are indicated by dotted and dashed lines, respectively. 

Numerical results, obtained from simulations at i?^, 
are shown in Fig. ^ After a short transient the system 
enters an intermediate regime, which extends up to sev- 
eral hundred time steps. Here the active sites form a 
single cluster and we observe power-law behavior with 
CDP exponents (dotted lines in Fig. This intermedi- 
ate regime is followed by a long crossover from CDP to 
DP, extending over almost two decades up to more than 
10* time steps, after which the system enters an asymp- 
totic DP regime (indicated by dashed lines in Fig. ^). 

Compared with ordinary DP lattice models, this cross- 
over regime is extremely long. We observed that by in- 
creasing / the crossover time can be reduced by more 
than one decade. Hence, for an experimental verification 
of DP, systems with high dissipation are more appropri- 
ate. The present experiments correspond to about 3000 
time steps (rows of beads); increasing this to about lO'' 
by using a longer inclined plane and smaller beads should 
yield DP behavior, provided that deviations of the exper- 
iment from the model do not increase with system size. 

The crossover from CDP to DP is illustrated in Fig. ||. 
Two avalanches are plotted on different scales. The left 
one represents a typical avalanche within the first few 
thousand time steps. As can be seen, the cluster appears 
to be compact on a lateral scale up to 100 lattice sites. 
However, as can be seen in the right panel of Fig. ||, 
after very long time the cluster breaks up into several 
branches, displaying the typical patterns of critical DP 
clusters. Thus, before measuring critical exponents, this 
feature has to be tested experimentally. To this end the 
DD experiment should be performed repeatedly at the 
critical tilt ip — ipQ. In most cases the avalanches will be 
small and compact. However, sometimes large avalanches 
will be generated which reach the bottom of the plate. 
If these avalanches are non-compact, we expect DP-type 
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FIG. 5. Typical clusters generated at criticality on small 
and large scales, illustrating the crossover from CDP to DP. 

asymptotic critical behavior. Only then is it worthwhile 
to optimize the experimental setup and to measure the 
critical exponents quantitatively. 
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